── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Registered S3 method overwritten by 'seriation':
method from
reorder.hclust vegan
Tip: ggsave to save the day your plots? 🦸📊💾
Viewing your plots in the RStudio Plots window is okay for practice
But don’t copy paste them to save them, as there is a much better way!
# Assign your plot to an R objectcontrol_barplot<-ps%>%ps_filter(case_control=="Control")%>%comp_barplot("Family", n_taxa =10, merge_other =FALSE)+coord_flip()+ggtitle("Control participants")# create a folder for figuresdir.create(here("practical-2/figs"))# Write the plot to a file, with ggsaveggsave( plot =control_barplot, filename =here("practical-2/figs/control_bars.png"), width =6.5, height =6.5, units ="in", dpi =300)
Be sure to carefully adjust the sizing and resolution of your plots for your paper or presentation!
Tip: Sizing plots in Quarto reports 📏📊📄
Again, viewing your plots in the RStudio Plots window is only okay for practice.
If you’re using Quarto reports, you can specify figure size and resolution in each chunk
e.g.
TODO display example chunk options examples directly in instructions - look up how to do this in quarto website
5.1 Organising your bar charts
Let’s look at all of the data from all diagnosis groups.
We can add the ggplot2 facet_wrap() to our plot, to separate the groups.
# A tibble: 2 × 6
term df sumsq meansq statistic p.value
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 diagnosis 2 1837. 919. 7.83 0.000745
2 Residuals 87 10203. 117. NA NA
TODO add posthoc pairwise comparisons
# TODO posthoc
Combine stats and plots in one go with ggstatsplot betweenstats
# TODO
6.2 Diversity
A true measure of ecosystem diversity (e.g. Shannon index) will consider the richness and evenness of the ecosystem.
Remember: a rich ecosystem predominated by one taxon is still intuitively a less diverse ecosystem than one with an even distribution of the same number of taxa.
We will use the effective Shannon diversity of genera \(e^H\) - because of its more intuitive interpretation.
ps<-ps_calc_diversity(ps, rank ="Genus", index ="shannon", varname ="Shannon_Genus")ps<-ps%>%ps_mutate(Effective_Shannon_Genus =exp(Shannon_Genus))
microViz provides a Shiny app ord_explore to interactively create and explore PCoA plots and other ordinations. See the code below to get started. But read the instructions first.
Instructions: a few things to try out
Colour the samples using the variables in the sample data
Select a few samples to view their composition on bar charts!
Change some ordination options:
Different rank of taxonomic aggregation
Different distances we’ve discussed
Copy the automatically generated code
Exit the app (press escape or click red 🛑 button in R console!)
Paste and run the code to recreate the ordination plot
Customise the plot: change colour scheme, title, etc.
Launch the app again with a different subset of the data
Practice using ps_filter etc.
e.g. plot the data of only the UC patients’ gut microbiota!
compute and colour points by an alpha diversity measure?
Beware: some important notes on interactive analysis
Unblock pop-ups: To allow the interactive analysis window to open in your browser, you may need to unblock pop-ups for your AMI instance address (check for messages about this after running the ord_explore command)
Slow UniFrac: UniFrac distances can be quite slow (over a minute) to calculate!
Filter to fewer samples and fewer taxa to speed it up (Before launching the app)
Many other distances: There are many distances available, feel free to try out ones we haven’t talked about
BUT:
You shouldn’t use a distance that you don’t understand in your actual work, even if the plot looks nice! 😉
A few of the distances might not work…
They are mostly implemented in the package vegan and I haven’t tested them all
Errors will appear in the RStudio R console
You can report to me any distances that don’t work (if you’re feeling helpful! 😇)
Other ordination methods: There are other ordination methods available in ord_explore
Try out PCA, principal components analysis, which does NOT use distances
We will not cover constrained and conditioned ordinations
If you are interested in e.g. RDA, you can look this up later
# fire up the shiny app# run these lines in your console (don't keep in script/notebook)ps%>%tax_filter(min_prevalence =2, verbose =FALSE)%>%# calculate new sample variables with dominant taxon (optional)ps_calc_dominant(rank ="Genus", none ="Mixed", other ="Other")%>%# launch a Shiny app in your web browser!ord_explore()
# Again, with different options# Run these lines in your consoleps%>%ps_filter(diagnosis=="UC")%>%tax_filter(min_prevalence =2, verbose =FALSE)%>%# calculate a few sample variables for interest (optional)ps_calc_dominant(rank ="Genus", none ="Mixed", other ="Other")%>%ps_calc_diversity(rank ="Genus", index ="shannon")%>%ps_calc_richness(rank ="Genus", index ="observed")%>%ord_explore()
10 PERMANOVA
Permutational multivariate analysis of variance.
ANOVA - analysis of variance (statistical modelling approach)
Multivariate - more than one dependent variable (multiple taxa!)
Permutational - statistical significance estimates obtained by shuffling the data many times
2024-06-28 10:47:48.227468 - Starting PERMANOVA with 99 perms with 1 processes
2024-06-28 10:47:48.245568 - Finished PERMANOVA
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 99
vegan::adonis2(formula = formula, data = metadata, permutations = n_perms, by = by, parallel = parall)
Df SumOfSqs R2 F Pr(>F)
diagnosis 2 1.5945 0.05962 2.7579 0.01 **
Residual 87 25.1499 0.94038
Total 89 26.7444 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Use more permutations for a more reliable p.value in your real work (slower)# Set a random seed number for reproducibility of this stochastic method
You can see from the model output that the p value, Pr(>F) is below 0.05. So there is good statistical evidence that the bacterial gut microbiota composition of C-section delivered infants has a different composition than vaginally delivered infants at 4 days of age.
10.1 Reporting PCoA and PERMANOVA methods
Your methodological choices matter, you should report what you did:
any relevant rare taxon filtering thresholds
the taxonomic rank of aggregation
the dissimilarity measure used to compute the pairwise distances
It’s probably a good idea to decide on a couple of appropriate distance measures up front for these tests, and report both (at least in supplementary material), as the choice of distance measure can affect results and conclusions!
10.2 Covariate-adjusted PERMANOVA
You can also adjust for covariates in PERMANOVA, and often should, depending on your study design. Let’s fit a more complex model, adjusting for sex and age.
The sequencing data gives us relative abundances, not absolute abundances.
The total number of reads sequenced per sample is an arbitrary total.
This leads to two main types of problem:
Interpretation caveats: see differential abundance section later
Statistical issues: taxon abundances are not independent, and may appear negatively correlated
These issues are worse with simpler ecosystems
Example: If one taxon blooms, the relative abundance of all other taxa will appear to decrease, even if they did not.
There is the same problem in theory with RNAseq data, but I suspect it is less bothersome because there are many more competing “species” of RNA transcript than there are bacterial species in even a very complex microbiome.The centered-log-ratio transformation (along with some other similar ratio transformations) are claimed to help with the statistical issues by transforming the abundances from the simplex to the real space.
TL;DR - the CLR transformation is useful for compositional microbiome data.
Practically, the CLR transformation involves finding the geometric mean of each sample
Then dividing abundance of each taxon in that sample by this geometric mean
Finally, you take the natural log of these ratios
ps%>%tax_sort(by =sum, at ="Genus", trans ="compositional", tree_warn =FALSE)%>%tax_agg(rank ="Genus")%>%tax_transform(trans ="clr", zero_replace ="halfmin", chain =TRUE)%>%comp_heatmap( samples =1:20, taxa =1:20, grid_lwd =2, name ="CLR", colors =heat_palette(sym =TRUE), tax_seriation ="Identity")
Interestingly, samples on the right of the plot (which tend to be UC patients) seem to have relatively more Escherichia/Shigella, and less Blautia, Faecalibacterium and Roseburia.
Wait, how to interpret these taxa loadings?
Cautiously
There are caveats and nuance to the interpretation of these plots, which are called PCA bi-plots
The relative length and direction of an arrow indicates how much that taxon contributes to the variation on each visible PC axis, e.g. Variation in Faecalibacterium abundance contributes quite a lot to variation along the PC1 axis.
The direction allows you to infer that samples positioned towards the left of the plot will tend to have higher relative abundance of Faecalibacterium than samples on the right of the plot.
Bacteroides variation contributes to both PC1 and PC2, as indicate by its high (negative) values on both axes.
Fancy circular bar charts?
We can make another kind of bar plot, using the PCA information to order our samples in a circular layout.
We probably want to filter out rare taxa, before performing some kinds of analysis.
Why remove rare taxa?
Rare taxa might sometimes be:
Sequencing errors
Statistically problematic
Biologically irrelevant
How to remove rare taxa?
What is rare? Two main concepts.
Low prevalence - taxon only detected in a small number of samples in your dataset.
Low abundance - relatively few reads assigned to that taxon (on average or in total)
Considering the impact of issues 1, 2, and 3, let’s say we are not interested in Species that occur in fewer than 2% of samples, and they have to have at least 10,000 reads in total across all samples.
ps%>%tax_filter(min_prevalence =2, min_total_abundance =100)%>%ntaxa()# after filtering
[1] 299
Wow so that would remove most of our unique taxa!
What is going on? Let’s make some plots!
# make table of summary statistics for the unique taxa in shao19psTaxaStats<-tibble( taxon =taxa_names(ps), prevalence =microbiome::prevalence(ps), total_abundance =taxa_sums(ps))
Warning: The `trans` argument of `sec_axis()` is deprecated as of ggplot2 3.5.0.
ℹ Please use the `transform` argument instead.
Some ggplot2 code
p
So most taxa have a low prevalence, and handful have way more reads than most.
Let’s label those points to check which taxa are the big time players.
p+ggrepel::geom_text_repel( data =function(df)filter(df, total_abundance>1500|prevalence>0.6), mapping =aes(label =taxon), size =2.5, min.segment.length =0, force =20, nudge_y =0.05)
Those taxa make sense for this dataset of gut microbiota samples.
Now let’s zoom in on the less abundant taxa by log-transforming the axes. We’ll also add lines indicating the thresholds of 2% prevalence and 1000 reads abundance.
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `Phylum = forcats::fct_explicit_na(Phylum, na_level = "Other")`.
Caused by warning:
! `fct_explicit_na()` was deprecated in forcats 1.0.0.
ℹ Please use `fct_na_value_to_level()` instead.
From the PCA loadings and barplots below, we have some strong suspicions about which taxa have a higher relative abundance in vaginally delivered infants than in c-section delivered infants, and vice versa, but we can also statistically test this. This is often called “differential abundance” (DA) testing, in the style of “differential expression” (DE) testing from the transcriptomics field.
This method simple method is borrowed from MaAsLin2
Note: they call the compositional transformation “Total Sum Scaling (TSS)”)
This is quite a straightforward method, so we will stick with this for today
But, many statistical methods have been developed for differential abundance analyses
Microbiome abundance data are quite awkward, statistically speaking, due to their sparseness and compositionality. Each successive method claims to handle some aspect of this awkwardness “better” than previous methods.
The aim is to have a method with adequate power to detect true associations, whilst controlling the type 1 error rate, the “false positive” reporting of associations that are not “truly” present.
Filter out the noise & interpret results with caution! use multiple testing corrections
Remember it’s all relative (abundance)
Try 2 or 3 methods and/or use same method as a previous study if replicating
Avoid Lefse and edgeR?
Beware: Not all methods allow covariate adjustment & few allow random effects (for time-series)
13.4 Now model all the taxa!?
We’re not normally interested in just one taxon!
It’s also hard to decide which taxonomic rank we are most interested in modelling!
Lower ranks like species or ASVs give better resolution but also more sparsity and classification uncertainty…
Higher ranks e.g. classes, could also be more powerful if you think most taxa within that class will follow a similar pattern.
So now we will fit a similar model for almost every taxon* at every rank from phylum to genus
*We’ll filter out species with a prevalence of less than 10%
# The code for `taxatree_models` is quite similar to tax_model.# However, you might need to run `tax_prepend_ranks` to ensure that each taxon at each rank is always unique.psModels<-ps%>%tax_prepend_ranks()%>%tax_transform("compositional", rank ="Genus", keep_counts =TRUE)%>%tax_filter(min_prevalence =0.1, undetected =0, use_counts =TRUE)%>%taxatree_models( type =lm, trans ="log2", trans_args =list(zero_replace ="halfmin"), ranks =c("Phylum", "Class", "Order", "Family", "Genus"), variables =c("IBD", "Female", "Age_Z"))
psModels
psExtra object - a phyloseq object with extra slots:
phyloseq-class experiment-level object
otu_table() OTU Table: [ 65 taxa and 90 samples ]
sample_data() Sample Data: [ 90 samples by 23 sample variables ]
tax_table() Taxonomy Table: [ 65 taxa by 6 taxonomic ranks ]
otu_get(counts = TRUE) [ 65 taxa and 90 samples ]
psExtra info:
tax_agg = "Genus" tax_trans = "compositional"
taxatree_models list:
Ranks: Phylum/Class/Order/Family/Genus
Why filter the taxa? It’s less likely that we are interested in rare taxa, and models of rare taxon abundances are more likely to be unreliable. Reducing the the number of taxa modelled also makes the process faster and makes visualizing the results easier!
Getting stats from the models
Next we will get a data.frame containing the regression coefficient estimates, test statistics and corresponding p values from all these regression models.
We have performed a lot of statistical tests here!
It is likely that we could find some significant p-values by chance alone.
We should correct for multiple testing / control the false discovery rate or family-wise error rate.
Instead of applying these adjustment methods across all taxa models at all ranks, the default behaviour is to control the family-wise error rate per taxonomic rank.
But how do we know which taxa are which nodes? We can create a labelled grey tree with taxatree_plotkey(). This labels only some of the taxa based on certain conditions that we specify.
set.seed(123)# label positionkey<-psStats%>%taxatree_plotkey( taxon_renamer =function(x)stringr::str_remove(x, "[PFG]: "),# conditions below, for filtering taxa to be labelledrank=="Phylum"|rank=="Genus"&prevalence>0.2# all phyla are labelled, and all genera with a prevalence of over 0.2)key
You can do more with these trees to customise them to your liking. See an extended tutorial here on the microViz website: including how to directly label taxa on the colored plots, change the layout and style of the trees, and even how to use a different regression modelling approach.
---title: "Practical 2 - microbiota analyses walkthrough"subtitle: "2024 NUTRIM microbiome & metabolome workshop"author: David Barnettdate: last-modifiedformat: htmlkeep-md: falsetheme: light: flatly dark: darklycss: ../../.css/instructions.cssembed-resources: truecode-block-border-left: truecode-block-bg: truetoc: truetoc-location: righttoc-depth: 4toc-expand: 1number-sections: truenumber-depth: 3fig-align: centerfig-dpi: 300fig-width: 7.5fig-height: 5fig-responsive: truecode-tools: truecode-fold: falsecode-link: truelightbox: autocache: false---## Intro**WORK IN PROGRESS** - this tutorial is currently just roughly adapted from previous content with a different datasetMore refinement and restructuring is planned.## Research questions 💪1. **Primary aim:** Does the bacterial gut microbiota composition of IBD-diagnosed patients differ from the control patients? a. **Diversity:** Is richness or diversity associated with IBD diagnosis? b. **Composition:** Does overall bacterial microbiota composition associated with IBD diagnosis? c. **Taxa:** Is the relative abundance of specific bacterial taxa (e.g. genera) associated with IBD diagnosis?2. **Activity:** Is current activity level associated with microbiota diversity, composition, or the relative abundance of specific taxa?3. **Medication:** Are IBD-related medications associated with microbiota diversity, composition, or the relative abundance of specific taxa?## Load R packages 📦```{r}library(here)library(tidyverse)library(broom)library(phyloseq)library(microViz)```## Read and clean dataRead the phyloseq we created in part 1.```{r}ps <-read_rds(file =here("data/papa2012/processed/papa12_phyloseq.rds"))```------------------------------------------------------------------------## Plot compositions- Stacked compositional bar charts using the `microViz` package, `comp_barplot()` function.- Let's start with a smaller subset of the data, just the control group.::: callout-tip### Tip: filtering phyloseq samplesWe can filter the samples like this, using the sample_data information```{r}ps %>%ps_filter(case_control =="Control") # similar to a dplyr filter!```If you forget a variable's levels, check with `table()` or `unique()````{r}ps@sam_data$case_control %>%table()```:::```{r, fig.width=5, fig.height=7, fig.align='center', out.width="70%"}ps %>% ps_filter(case_control == "Control") %>% comp_barplot("Family", n_taxa = 10, merge_other = FALSE) + coord_flip() + ggtitle("Control participants")```::: {.callout-tip collapse="true"}### Tip: `ggsave` to save ~~the day~~ your plots? 🦸📊💾- Viewing your plots in the RStudio Plots window is okay for practice- But don't copy paste them to save them, as there is a much better way!```{r}# Assign your plot to an R objectcontrol_barplot <- ps %>%ps_filter(case_control =="Control") %>%comp_barplot("Family", n_taxa =10, merge_other =FALSE) +coord_flip() +ggtitle("Control participants")# create a folder for figuresdir.create(here("practical-2/figs"))# Write the plot to a file, with ggsaveggsave(plot = control_barplot, filename =here("practical-2/figs/control_bars.png"),width =6.5, height =6.5, units ="in", dpi =300)``````{r, include=FALSE}unlink(here("practical-2/figs"), recursive = TRUE)```**Be sure to carefully adjust the sizing and resolution of your plots for your paper or presentation!**:::::: {.callout-tip collapse="true"}### Tip: Sizing plots in Quarto reports 📏📊📄- Again, viewing your plots in the RStudio Plots window is only okay for practice.- If you're using Quarto reports, you can specify figure size and resolution in each chunke.g. TODO display example chunk options examples directly in instructions - look up how to do this in quarto website:::### Organising your bar charts- Let's look at all of the data from all diagnosis groups.- We can add the ggplot2 `facet_wrap()` to our plot, to separate the groups.```{r}#| fig-width: 7#| fig-height: 5ps %>%comp_barplot("Family", n_taxa =10, merge_other =FALSE) +facet_wrap(facets =vars(diagnosis), scales ="free") +coord_flip() ```::: {.callout-note collapse="true"}### More barchart resources:More examples/tutorial of visualizing microbiota compositions using stacked bar charts can be found here: <https://david-barnett.github.io/microViz/articles/web-only/compositions.html>:::------------------------------------------------------------------------## ⚖️ Ecosystem diversity- How rich and diverse is the bacterial microbiota of each sample? - Does this richness or diversity differ by diagnosis group?### Richness- The simplest richness measure is just counting, a.k.a. "Observed Richness".- Let's compute the observed richness of genera.- Note: `ps_calc_richness()` computes the index for each sample and adds it to your sample_data```{r}ps <-ps_calc_richness(ps, rank ="Genus", index ="observed", varname ="N_genera")```::: panel-tabset#### Plot```{r}ps %>%samdat_tbl() %>%ggplot(aes(x = N_genera, y = diagnosis, color = diagnosis)) +geom_boxplot(outliers =FALSE) +geom_jitter(height =0.2) +theme_classic()```#### Linear regression / ANOVATODO decide what to show exactly - posthoc tests after anova? or lm coefficients```{r}richness_lm <-lm(data =samdat_tbl(ps), formula = N_genera ~ diagnosis)anova(richness_lm)```Get a tidy tibble of results with the `tidy()` function from `broom` package.```{r}richness_stats <- broom::tidy(anova(richness_lm))richness_stats```TODO add posthoc pairwise comparisons```{r}# TODO posthoc```#### `ggstatsplot`Combine stats and plots in one go with ggstatsplot betweenstats```{r}# TODO```:::### Diversity- A true measure of ecosystem diversity (e.g. Shannon index) will consider the richness and evenness of the ecosystem.- Remember: a rich ecosystem predominated by one taxon is still intuitively a less diverse ecosystem than one with an even distribution of the same number of taxa.- We will use the effective Shannon diversity of genera $e^H$ - because of its more intuitive interpretation.```{r}ps <-ps_calc_diversity(ps, rank ="Genus", index ="shannon", varname ="Shannon_Genus")ps <- ps %>%ps_mutate(Effective_Shannon_Genus =exp(Shannon_Genus))```::: panel-tabset#### Plot```{r}ps %>%samdat_tbl() %>%ggplot(aes(x = Effective_Shannon_Genus, y = diagnosis, color = diagnosis)) +geom_boxplot(outliers =FALSE) +geom_jitter(height =0.2) +theme_classic()```#### Linear Regression / ANOVA```{r}eShannon_lm <-lm(data =samdat_tbl(ps), formula = Effective_Shannon_Genus ~ diagnosis)anova(eShannon_lm)summary(eShannon_lm)```Get a tidy tibble of results with the `tidy()` function from `broom` package.```{r}eShannon_stats <- broom::tidy(anova(eShannon_lm))eShannon_stats```TODO add posthoc pairwise comparisons```{r}# TODO posthoc```#### `ggstatsplot`Combine stats and plots in one go with ggstatsplot betweenstats```{r}# TODO```#### Covariate adjusted```{r}```:::```{r}###```## Dissimilarity ### Dissimilarity measures- Calculating the dissimilarities between pairs of microbiome samples- We talked about these commonly-used dissimilarity measures in the lecture. - **Binary Jaccard** - presence-absence - **Bray-Curtis** - abundance-weighted - **UniFrac** distances (unweighted / weighted / generalised)```{r}ps %>%tax_agg(rank ="Genus") %>%tax_transform("binary") %>%# converts counts to absence/presence: 0/1dist_calc(dist ="jaccard")```::: {.callout-important collapse="true"}### "Binary" Jaccard `010101`- Remember to run a "binary" transform on your data before computing "jaccard" distance.- There is a quantitative form of the Jaccard distance, which is the default behaviour!- But the qualitative (presence/absence) version is mostly used in microbial ecology.- If you want an abundance-weighted ecological dissimilarity, use Bray-Curtis!:::- We now have our pairwise dissimilarities! 🎉- A distance matrix is attached as an extra part on the original phyloseq object::: {.callout-note collapse="true"}### Dissimilarity or distance?- These terms are often used interchangeably- You will find dissimilarities in a distance matrix- But if you want to be pedantic a true "distance metric" d, must satisfy 3 properties: 1. Identity of indiscernibles: For any samples $x$ and $y$, $d(x, y) = 0$ if and only if $x = y$ 2. Symmetry: For any samples $x$ and $y$, $d(x, y) = d(y, x)$ 3. Triangle inequality: For any samples $x$, $y$, and $z$, $d(x, z) ≤ d(x, y) + d(y, z)$- (3) can be interpreted as: "the direct path between two points must be at least as short as any detour"- This is not true for e.g. Bray-Curtis, but in practice it is very rarely problematic.:::- The object is now class `psExtra` (created by microViz)- A psExtra also stores info about the aggregation and transformations you performed```{r}distances <- ps %>%tax_agg(rank ="Genus") %>%tax_transform("binary") %>%dist_calc(dist ="jaccard") %>%dist_get()```You can extract the distance matrix with dist_get.```{r}as.matrix(distances)[1:5, 1:5]```The Binary Jaccard dissimilarities range between 0 (identical) and 1 (no shared genera).```{r}range(as.matrix(distances))```------------------------------------------------------------------------## Ordination- What can we do with these dissimilarities? 🤔- We can make an ordination! 💡- Conceptually, ordination refers to a process of ordering things (in our case: samples).- Similar samples are placed closer to each other, and dissimilar samples are placed further away.### PCoAPrincipal Co-ordinates Analysis is one kind of ordination.- PCoA takes a distance matrix and finds new dimensions (a co-ordinate system, if you like)- The new dimensions are created with the aim to preserve the original distances between samples- And to capture the majority of this distance information in the first dimensions- This makes it easier to visualize the patterns in your data, in 2D or 3D plots 👀::: {.callout-note collapse="true"}### For more info, see "GUSTAME"- There is helpful info about ordination methods including PCoA on the GUide to STatistical Analysis in Microbial Ecology website (GUSTA ME). <https://sites.google.com/site/mb3gustame/dissimilarity-based-methods/principal-coordinates-analysis>- This website covers a lot of other topics too, which may be interesting for you to read at a later date if you'll work on microbiome analysis.:::```{r}ps %>%tax_transform(trans ="identity", rank ="Genus") %>%dist_calc(dist ="bray") %>%ord_calc(method ="PCoA") %>%ord_plot(alpha =0.6, size =2) +theme_classic(12) +coord_fixed(0.7)```To get a little insight into what has happened here, we can colour each sample according to its dominant (most abundant) genus.```{r}ps %>%ps_calc_dominant(rank ="Genus", none ="Mixed", other ="Other", var ="Dominant Genus") %>%tax_transform(trans ="identity", rank ="Genus") %>%dist_calc(dist ="bray") %>%ord_calc(method ="PCoA") %>%ord_plot(color ="Dominant Genus", alpha =0.6, size =2) +scale_color_brewer(palette ="Dark2") +theme_classic(12) +coord_fixed(0.7)```## Interactive ordination!`microViz` provides a Shiny app `ord_explore` to interactively create and explore PCoA plots and other ordinations.See the code below to get started.But read the instructions first.::: {.callout-tip collapse="true"}### **Instructions:** a few things to try out- Colour the samples using the variables in the sample data- Select a few samples to view their composition on bar charts!- Change some ordination options: - Different rank of taxonomic aggregation - Different distances we've discussed- Copy the automatically generated code - Exit the app (press escape or click red 🛑 button in R console!) - Paste and run the code to recreate the ordination plot - Customise the plot: change colour scheme, title, etc.- Launch the app again with a different subset of the data - Practice using `ps_filter` etc. - e.g. plot the data of only the UC patients' gut microbiota! - compute and colour points by an alpha diversity measure?:::::: {.callout-warning collapse="true"}### **Beware: some important notes on interactive analysis**- **Unblock pop-ups**: To allow the interactive analysis window to open in your browser, you may need to unblock pop-ups for your AMI instance address (check for messages about this after running the ord_explore command)- **Slow UniFrac**: UniFrac distances can be quite slow (over a minute) to calculate! - Filter to fewer samples and fewer taxa to speed it up (Before launching the app)- **Many other distances**: There are many distances available, feel free to try out ones we haven't talked about - **BUT**: - You shouldn't use a distance that you don't understand in your actual work, even if the plot looks nice! 😉 - A few of the distances might not work... - They are mostly implemented in the package `vegan` and I haven't tested them all - Errors will appear in the RStudio R console - You can report to me any distances that don't work (if you're feeling helpful! 😇)- **Other ordination methods**: There are other ordination methods available in `ord_explore` - Try out PCA, principal **components** analysis, which does NOT use distances - We will not cover constrained and conditioned ordinations - If you are interested in e.g. RDA, you can look this up later - See the [Guide to Statistical Analysis in Microbial Ecology](https://sites.google.com/site/mb3gustame/constrained-analyses/redundancy-analysis):::``` r# fire up the shiny app# run these lines in your console (don't keep in script/notebook)ps %>%tax_filter(min_prevalence =2, verbose =FALSE) %>%# calculate new sample variables with dominant taxon (optional)ps_calc_dominant(rank ="Genus", none ="Mixed", other ="Other") %>%# launch a Shiny app in your web browser!ord_explore()`````` r# Again, with different options# Run these lines in your consoleps %>%ps_filter(diagnosis =="UC") %>%tax_filter(min_prevalence =2, verbose =FALSE) %>%# calculate a few sample variables for interest (optional)ps_calc_dominant(rank ="Genus", none ="Mixed", other ="Other") %>%ps_calc_diversity(rank ="Genus", index ="shannon") %>%ps_calc_richness(rank ="Genus", index ="observed") %>%ord_explore()```------------------------------------------------------------------------## PERMANOVAPermutational multivariate analysis of variance.- ANOVA - analysis of variance (statistical modelling approach)- Multivariate - more than one dependent variable (multiple taxa!)- Permutational - statistical significance estimates obtained by shuffling the data many times::: {.callout-note collapse="true"}### For more details on PERMANOVA- See this excellent book chapter by Marti Anderson on PERMANOVA: <https://onlinelibrary.wiley.com/doi/full/10.1002/9781118445112.stat07841>- Sometimes PERMANOVA is called NP-MANOVA (non-parametric MANOVA)- e.g. on the GUide to STatistical Analysis in Microbial Ecology [website](https://sites.google.com/site/mb3gustame/hypothesis-tests/manova/npmanova).:::- **TLDR:** Are those groups on the PCoA actually different??```{r}ps %>%tax_filter(min_prevalence =2, verbose =FALSE) %>%tax_agg(rank ="Genus") %>%dist_calc(dist ="bray") %>%ord_calc(method ="PCoA") %>%ord_plot(alpha =0.6, size =2, color ="diagnosis") +theme_classic(12) +coord_fixed(0.7) +stat_ellipse(aes(color = diagnosis)) +scale_color_brewer(palette ="Set1")``````{r}ps %>%tax_filter(min_prevalence =2, verbose =FALSE) %>%tax_agg(rank ="Genus") %>%dist_calc(dist ="bray") %>%dist_permanova(variables ="diagnosis", n_perms =99, seed =123) %>%perm_get()# Use more permutations for a more reliable p.value in your real work (slower)# Set a random seed number for reproducibility of this stochastic method```You can see from the model output that the p value, Pr(\>F) is below 0.05. So there is good statistical evidence that the bacterial gut microbiota composition of C-section delivered infants has a different composition than vaginally delivered infants at 4 days of age.### Reporting PCoA and PERMANOVA methods- Your methodological choices matter, you should report what you did: - any relevant rare taxon filtering thresholds - the taxonomic rank of aggregation - the dissimilarity measure used to compute the pairwise distancesIt's probably a good idea to decide on a couple of appropriate distance measures up front for these tests, and report both (at least in supplementary material), as the choice of distance measure can affect results and conclusions!### Covariate-adjusted PERMANOVAYou can also adjust for covariates in PERMANOVA, and often should, depending on your study design.Let's fit a more complex model, adjusting for sex and age.```{r}ps %>%tax_filter(min_prevalence =2, verbose =FALSE) %>%tax_agg(rank ="Genus") %>%dist_calc(dist ="bray") %>%dist_permanova(variables =c("diagnosis", "gender", "age_years"),n_perms =99, seed =111 ) %>%perm_get()# Use more permutations for a more reliable p.value in your real work (slower)# Set a random seed number for reproducibility of this stochastic method```## PCA- Principal **Components** Analysis.- For practical purposes, PCA is quite similar to Principal Co-ordinates Analysis.- In fact, PCA produces equivalent results to PCoA with Euclidean distances.::: {.callout-tip collapse="true"}### "Help, what are Euclidean distances?"- Euclidean distances are essentially a generalization of Pythagoras' theorem to more dimensions.- In our data every taxon is a feature, a dimension, on which we calculate Euclidean distances.**Pythagoras' theorem:**$$c = \sqrt{a^2 + b^2}$$**Euclidean distance:**$$d\left(p, q\right) = \sqrt{\sum _{i=1}^{n_{taxa}} \left( p_{i}-q_{i}\right)^2 }$$- Distance $d$ between samples $p$ and $q$, with $n$ taxa.:::::: panel-tabset#### Euclidean PCoA```{r}#| code-fold: trueps %>%tax_agg(rank ="Genus") %>%dist_calc(dist ="euclidean") %>%ord_calc(method ="PCoA") %>%ord_plot(alpha =0.6, size =2) +geom_rug(alpha =0.1)```#### PCA on counts```{r}#| code-fold: trueps %>%tax_agg(rank ="Genus") %>%ord_calc(method ="PCA") %>%ord_plot(alpha =0.6, size =2) +geom_rug(alpha =0.1)```:::**Problems with PCA (or PCoA with Euclidean Distances) on microbiota data**- These plots look weird! most samples bunch in the middle, with spindly projections..- Sensitive to sparsity (double-zero problem) --\> filter rare taxa- Excessive emphasis on high-abundance taxa --\> log transform features first------------------------------------------------------------------------## Log transformations, and CLR- First let's look at the abundance again, this time with heatmaps.- Each column is a sample (from a child), and each row is a taxon.```{r, fig.height=4, fig.width=6}ps %>% tax_sort(by = sum, at = "Genus", trans = "compositional", tree_warn = FALSE) %>% tax_transform(trans = "compositional", rank = "Genus") %>% comp_heatmap(samples = 1:20, taxa = 1:20, name = "Proportions", tax_seriation = "Identity")```- Even though we have picked the top 20 most abundant genera, there are still a lot of zeros- **Problem:** We need to deal with the zeros, because `log(0)` is undefined.- **Solution:** add a small amount to every value (or just every zero), before applying the log transformation.- This small value is often called a pseudo-count.::: {.callout-tip collapse="true"}### What value should we use for the pseudo-count?- One easy option is to just add a count of 1- Another popular option is to add half of the smallest observed real value (from across the whole dataset)- In general, for zero replacement, keep it simple and **record your approach**:::### Centered Log Ratio transformation:**Remember**, [Microbiome Datasets Are Compositional: And This Is Not Optional.](https://doi.org/10.3389/fmicb.2017.02224)::: {.callout-note collapse="true"}### More details on the "CoDa" problem:- The sequencing data gives us relative abundances, not absolute abundances.- The total number of reads sequenced per sample is an arbitrary total.**This leads to two main types of problem:**- Interpretation caveats: see differential abundance section later- Statistical issues: taxon abundances are not independent, and may appear negatively correlated- These issues are worse with simpler ecosystemsExample: If one taxon blooms, the relative abundance of all other taxa will appear to decrease, even if they did not.*There is the same problem in theory with RNAseq data, but I suspect it is less bothersome because there are many more competing "species" of RNA transcript than there are bacterial species in even a very complex microbiome.* *The centered-log-ratio transformation (along with some other similar ratio transformations) are claimed to help with the statistical issues by transforming the abundances from the simplex to the real space.*:::**TL;DR - the CLR transformation is useful for compositional microbiome data.**- Practically, the CLR transformation involves finding the geometric mean of each sample- Then dividing abundance of each taxon in that sample by this geometric mean- Finally, you take the natural log of these ratios```{r, fig.height=3, fig.width=6}ps %>% tax_sort(by = sum, at = "Genus", trans = "compositional", tree_warn = FALSE) %>% tax_agg(rank = "Genus") %>% tax_transform(trans = "clr", zero_replace = "halfmin", chain = TRUE) %>% comp_heatmap( samples = 1:20, taxa = 1:20, grid_lwd = 2, name = "CLR", colors = heat_palette(sym = TRUE), tax_seriation = "Identity" )```### PCA on CLR-transformed taxa```{r}ps %>%tax_filter(min_prevalence =2, verbose =FALSE) %>%tax_transform(rank ="Genus", trans ="clr", zero_replace ="halfmin") %>%ord_calc(method ="PCA") %>%ord_plot(alpha =0.6, size =2, color ="diagnosis") +theme_classic(12) +coord_fixed(0.7)```- After the CLR transformation, the plot looks better- We can see a pattern where the gut microbiomes of infants cluster by birth mode#### So why is PCA interesting for us?- Principal components are built directly from a (linear) combination of the original features.- That means we know how much each taxon contributes to each PC dimension- We can plot this information (loadings) as arrows, alongside the sample points```{r}pca <- ps %>%tax_filter(min_prevalence =2, verbose =FALSE) %>%tax_transform(rank ="Genus", trans ="clr", zero_replace ="halfmin") %>%ord_calc(method ="PCA") %>%ord_plot(alpha =0.6, size =2, color ="diagnosis", plot_taxa =1:6, tax_vec_length =0.5,tax_lab_style =tax_lab_style(type ="text", max_angle =90, aspect_ratio =1,size =3, fontface ="bold" ), ) +theme_classic(12) +coord_fixed(ratio =1, xlim =c(-3, 3), ylim =c(-3, 3), clip ="off")pca```Interestingly, samples on the right of the plot (which tend to be UC patients) seem to have relatively more *Escherichia/Shigella*, and less *Blautia*, *Faecalibacterium* and *Roseburia*.::: {.callout-important collapse="true"}### Wait, how to interpret these taxa loadings?##### Cautiously- There are caveats and nuance to the interpretation of these plots, which are called PCA bi-plots- You can read more here: <https://sites.google.com/site/mb3gustame/indirect-gradient-analysis/principal-components-analysis>##### In general:The relative length and direction of an arrow indicates how much that taxon contributes to the variation on each visible PC axis, e.g. Variation in *Faecalibacterium* abundance contributes quite a lot to variation along the PC1 axis.The direction allows you to infer that samples positioned towards the left of the plot will tend to have higher relative abundance of *Faecalibacterium* than samples on the right of the plot.*Bacteroides* variation contributes to both PC1 and PC2, as indicate by its high (negative) values on both axes.:::::: {.callout-note collapse="true"}### Fancy circular bar charts?We can make another kind of bar plot, using the PCA information to order our samples in a circular layout.```{r}iris <- ps %>%tax_filter(min_prevalence =2, verbose =FALSE) %>%tax_transform(rank ="Genus", trans ="clr", zero_replace ="halfmin") %>%ord_calc(method ="PCA") %>%ord_plot_iris(tax_level ="Genus", n_taxa =12, other ="Other",anno_colour ="diagnosis",anno_colour_style =list(alpha =0.6, size =0.6, show.legend =FALSE) )``````{r, fig.height=5, fig.width=10}patchwork::wrap_plots(pca, iris, nrow = 1, guides = "collect")```:::::: {.callout-note collapse="true"}### Notes on filtering rare taxaWe probably want to filter out **rare** taxa, before performing some kinds of analysis.##### Why remove rare taxa?**Rare taxa might sometimes be:**1. Sequencing errors2. Statistically problematic3. Biologically irrelevant##### How to remove rare taxa?**What is rare?** Two main concepts.- Low **prevalence** - taxon only detected in a small number of samples in your dataset.- Low **abundance** - relatively few reads assigned to that taxon (on average or in total)Considering the impact of issues 1, 2, and 3, let's say we are not interested in Species that occur in fewer than 2% of samples, and they have to have at least 10,000 reads in total across all samples.```{r}ntaxa(ps) # before filtering``````{r}ps %>%tax_filter(min_prevalence =2, min_total_abundance =100) %>%ntaxa() # after filtering```- Wow so that would remove **most** of our unique taxa!- What is going on? Let's make some plots!```{r}# make table of summary statistics for the unique taxa in shao19psTaxaStats <-tibble(taxon =taxa_names(ps),prevalence = microbiome::prevalence(ps),total_abundance =taxa_sums(ps))``````{r}#| code-fold: true#| code-summary: Some ggplot2 codep <- psTaxaStats %>%ggplot(aes(total_abundance, prevalence)) +geom_point(alpha =0.5) +geom_rug(alpha =0.1) +scale_x_continuous(labels = scales::label_number(), name ="Total Abundance") +scale_y_continuous(labels = scales::label_percent(), breaks = scales::breaks_pretty(n =9),name ="Prevalence (%)",sec.axis =sec_axis(trans =~ . *nsamples(ps), breaks = scales::breaks_pretty(n =9),name ="Prevalence (N samples)" ) ) +theme_bw()p```So most taxa have a low prevalence, and handful have way more reads than most.Let's label those points to check which taxa are the big time players.```{r}p + ggrepel::geom_text_repel(data =function(df) filter(df, total_abundance >1500| prevalence >0.6),mapping =aes(label = taxon), size =2.5, min.segment.length =0, force =20, nudge_y =0.05)```Those taxa make sense for this dataset of gut microbiota samples.Now let's zoom in on the less abundant taxa by log-transforming the axes.We'll also add lines indicating the thresholds of 2% prevalence and 1000 reads abundance.```{r}#| code-fold: true#| code-summary: Some more ggplot2 codepsTaxaStats %>%ggplot(aes(x = total_abundance, y = prevalence)) +geom_vline(xintercept =10, color ="red", linetype ="dotted") +geom_hline(yintercept =2/100, color ="red", linetype ="dotted") +geom_point(alpha =0.5) +geom_rug(alpha =0.1) +scale_x_log10(labels = scales::label_number(), name ="Total Abundance") +scale_y_log10(labels = scales::label_percent(), breaks = scales::breaks_log(n =9),name ="Prevalence (%)",sec.axis =sec_axis(trans =~ . *nsamples(ps), breaks = scales::breaks_log(n =9),name ="Prevalence (N samples)" ) ) +theme_bw()```- We can break this down by phylum if we add the taxonomic table information.```{r, fig.height = 5, fig.width=8}#| code-fold: true#| code-summary: A lot more ggplot2 code!# don't worry about this code if it's confusing, just focus on the plot outputps %>% tax_table() %>% as.data.frame() %>% as_tibble(rownames = "taxon") %>% left_join(psTaxaStats, by = "taxon") %>% add_count(Phylum, name = "phylum_count", sort = TRUE) %>% mutate(Phylum = factor(Phylum, levels = unique(Phylum))) %>% # to fix facet order mutate(Phylum = forcats::fct_lump_n(Phylum, n = 5)) %>% mutate(Phylum = forcats::fct_explicit_na(Phylum, na_level = "Other")) %>% ggplot(aes(total_abundance, prevalence)) + geom_vline(xintercept = 10, color = "red", linetype = "dotted") + geom_hline(yintercept = 2 / 100, color = "red", linetype = "dotted") + geom_point(alpha = 0.5, size = 1) + geom_rug(alpha = 0.2) + scale_x_log10( labels = scales::label_log(), breaks = scales::breaks_log(n = 5), name = "Total Abundance" ) + scale_y_log10( labels = scales::label_percent(), breaks = scales::breaks_log(n = 9), name = "Prevalence (%)", sec.axis = sec_axis( trans = ~ . * nsamples(shao19), breaks = scales::breaks_log(n = 9), name = "Prevalence (N samples)" ) ) + facet_wrap("Phylum") + theme_bw(10)```**How you pick a threshold, depends on what analysis method you are filtering for!**- alpha diversity: do not filter- beta diversity: relevance of threshold depends on your distance measure- differential abundance testing: stringent filtering, prevalence \>5%, \>10%? :::------------------------------------------------------------------------## Taxon statsFrom the PCA loadings and barplots below, we have some strong suspicions about which taxa have a higher relative abundance in vaginally delivered infants than in c-section delivered infants, and vice versa, but we can also statistically test this.This is often called "differential abundance" (DA) testing, in the style of "differential expression" (DE) testing from the transcriptomics field.```{r}ps %>%comp_barplot(tax_level ="Genus", n_taxa =12, facet_by ="diagnosis", label =NULL, ) +coord_flip() +theme(axis.ticks.y =element_blank())```### Model one taxon- We will start by creating a linear regression model for one genus, Bacteroides.- We will transform the count data by first making it proportions, and then taking a base 2 logarithm, log2, after adding a pseudocount.```{r}bacteroidesRegression1 <- ps %>%tax_transform("compositional", rank ="Genus") %>%tax_model(type ="lm", rank ="Genus",trans ="log2", trans_args =list(zero_replace ="halfmin"),taxa ="Bacteroides", variables ="case_control",return_psx =FALSE ) %>%pluck(1)```- Looking at the regression results```{r}summary(bacteroidesRegression1)``````{r}broom::tidy(bacteroidesRegression1, conf.int =TRUE)```### Covariate-adjusted modelWe can fit a model with covariates, as we did for PERMANOVA- We will convert the categorical variables into indicator (dummy) variables - We will scale the continuous covariates to 0 mean and SD 1 (z-scores)- You'll see this will make our subsequent plots easier to interpret later```{r}ps <- ps %>%ps_mutate(IBD =if_else(case_control =="Case", true =1, false =0),Female =if_else(gender =="female", true =1, false =0),Age_Z =scale(age_years, center =TRUE, scale =TRUE) )``````{r}bacteroidesRegression2 <- ps %>%tax_transform("compositional", rank ="Genus") %>%tax_model(type ="lm", rank ="Genus", taxa ="Bacteroides",trans ="log2", trans_args =list(zero_replace ="halfmin"),variables =c("IBD", "Female", "Age_Z"),return_psx =FALSE ) %>%pluck(1)```- Looking at the regression results```{r}summary(bacteroidesRegression2)broom::tidy(bacteroidesRegression2, conf.int =TRUE)```### There are many DA methods!- This method simple method is borrowed from MaAsLin2 - Note: they call the compositional transformation "Total Sum Scaling (TSS)")- This is quite a straightforward method, so we will stick with this for today- But, many statistical methods have been developed for differential abundance analysesMicrobiome abundance data are quite awkward, statistically speaking, due to their sparseness and compositionality.Each successive method claims to handle some aspect of this awkwardness "better" than previous methods.The aim is to have a method with adequate power to detect true associations, whilst controlling the type 1 error rate, the "false positive" reporting of associations that are not "truly" present.Results are surprisingly inconsistent across the different methods, as demonstrated this year in a [fascinating analysis by Jacob Nearing and colleagues](https://www.nature.com/articles/s41467-022-28034-z).#### So, what to do?- Filter out the noise & interpret results with caution! use multiple testing corrections- Remember it's all relative (abundance)- Try 2 or 3 methods and/or use same method as a previous study if replicating - Avoid Lefse and edgeR? - Beware: Not all methods allow covariate adjustment & few allow random effects (for time-series)### Now model all the taxa!?1. We're not normally interested in just one taxon!2. It's also hard to decide which taxonomic rank we are most interested in modelling! - Lower ranks like species or ASVs give better resolution but also more sparsity and classification uncertainty... - Higher ranks e.g. classes, could also be more powerful if you think most taxa within that class will follow a similar pattern.- So now we will fit a similar model for almost every taxon\* at every rank from phylum to genus- \*We'll filter out species with a prevalence of less than 10%```{r}#| warning: false# The code for `taxatree_models` is quite similar to tax_model.# However, you might need to run `tax_prepend_ranks` to ensure that each taxon at each rank is always unique.psModels <- ps %>%tax_prepend_ranks() %>%tax_transform("compositional", rank ="Genus", keep_counts =TRUE) %>%tax_filter(min_prevalence =0.1, undetected =0, use_counts =TRUE) %>%taxatree_models(type = lm,trans ="log2", trans_args =list(zero_replace ="halfmin"),ranks =c("Phylum", "Class", "Order", "Family", "Genus"),variables =c("IBD", "Female", "Age_Z") )``````{r}psModels```*Why filter the taxa? It's less likely that we are interested in rare taxa, and models of rare taxon abundances are more likely to be unreliable. Reducing the the number of taxa modelled also makes the process faster and makes visualizing the results easier!*#### Getting stats from the modelsNext we will get a data.frame containing the regression coefficient estimates, test statistics and corresponding p values from all these regression models.```{r}psStats <-taxatree_models2stats(psModels)psStats``````{r}psStats %>%taxatree_stats_get()```#### Adjusting p values- We have performed a lot of statistical tests here!- It is likely that we could find some significant p-values by chance alone.- We should correct for multiple testing / control the false discovery rate or family-wise error rate.*Instead of applying these adjustment methods across all taxa models at all ranks, the default behaviour is to control the family-wise error rate per taxonomic rank.*```{r}psStats <- psStats %>%taxatree_stats_p_adjust(method ="BH", grouping ="rank")# notice the new variablepsStats %>%taxatree_stats_get()```### Plot all the taxatree_stats!- `taxatree_plots()` allows you to plot statistics from all of the taxa models onto a tree layout (e.g. point estimates and significance).- The taxon model results are organised by rank, radiating out from the central root node - e.g. from Phyla around the center to Genus in the outermost ring.`taxatree_plots()` itself returns a list of plots, which you can arrange into one figure with the [`patchwork`](https://patchwork.data-imaginist.com/) package for example (and/or [`cowplot`](https://wilkelab.org/cowplot/articles/plot_grid.html)).```{r, fig.width=6, fig.height=5}psStats %>% taxatree_plots(node_size_range = c(1, 3), sig_stat = "p.adj.BH.rank") %>% patchwork::wrap_plots(ncol = 2, guides = "collect")```#### Taxatree KeyBut how do we know which taxa are which nodes?We can create a labelled grey tree with `taxatree_plotkey()`.This labels only some of the taxa based on certain conditions that we specify.```{r fig.height=4, fig.width=4.5, warning=FALSE}set.seed(123) # label positionkey <- psStats %>% taxatree_plotkey( taxon_renamer = function(x) stringr::str_remove(x, "[PFG]: "), # conditions below, for filtering taxa to be labelled rank == "Phylum" | rank == "Genus" & prevalence > 0.2 # all phyla are labelled, and all genera with a prevalence of over 0.2 )key```You can do more with these trees to customise them to your liking.See an extended tutorial [here on the microViz website](https://david-barnett.github.io/microViz/articles/web-only/modelling-taxa.html#plot-all-the-taxatree_stats): including how to directly label taxa on the colored plots, change the layout and style of the trees, and even how to use a different regression modelling approach.```{r}# try it out!```## Next! ⏩## Session info<details>```{r}sessioninfo::session_info()```</details>
